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ABSTRACT 


The use of Doubly-Fed Induction Generators (DFIG) for wind energy conversion is 
addressed in this thesis. It is well known that when the stator is connected to the electric 


grid, the rotor voltage can control both mechanical torque and reactive electric power. 


To guarantee efficient wind energy conversion, it is important to research and 
design more advanced control schemes. In this thesis, we first review the basic theory 
behind DFIGs and Adaptive Control. Next we design an adaptive controller for a wind 
turbine using a DFIG and model and simulate the system. In order to create a valid 
assessment on the results of this method, we compare the system’s performance with a 
standard control scheme based on proportional integral (PI) controllers as proposed in 


standard approaches. 
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EXECUTIVE SUMMARY 


The possibility of fossil fuel shortages in the near future and environmental concerns 
about power generation has lead to an increasing demand for renewable energy (RE) 
installations worldwide. Among them, wind energy (WE) is one of the fastest growing. 
This form of energy production is suitable for electricity production since wind is 
available everywhere. However, the random nature of the wind and the need to provide a 
stable source for the grid calls for the design of effective control systems for wind energy 
conversion systems (WECS). This would improve both the reliability of the energy 


captured and the overall efficiency of the system. 


Many large wind power generators today are based on the doubly-fed induction 
generator (DFIG). The DFIG is capable of controlling both the active and reactive power 
simultaneously and keeping the generator at its maximum efficiency, even in the 
presence of wind speed variations. To take advantage of these features, the DFIG requires 


proper control in terms of its power electronics and computer algorithms. 


Standard controllers, such as proportional integral (PI) controllers for pulse width 
modulation (PWM) converters at the rotor and source side of the system are the most 


commonly used to provide independent control of active and reactive power output. 


A different control approach based on adaptive control theory is investigated in 
this thesis. The motivation for this research is the fact that WECS operate in 
environments with extreme variations in the operating conditions, and their operation 
depends on the wind speed, which, in general is a random process that cannot be 
controlled. Adaptive control systems, on the other hand, have the ability to “adapt” to 
variations in the system’s parameters, which is a very attractive feature. Among the 
various kinds of adaptive control methods, the one examined in this thesis is direct model 
reference adaptive control (DMRAC). The basic idea behind this approach is to make the 
system behave as closely as possible to a reference model in the presence of uncertain 
and/or time-varying dynamics of the system to be controlled. To achieve this, the system 


is modeled in terms of a vector of parameters 6,(t) that may be known or approximately 


XV 


known. This parameter vector is updated by an adaptive law using input and output data. 


The goal of the controller is to continuously estimate the parameters creating a vector 
6.(t) to drive the error of the output vector z(t) between the reference model and the 


actual plant to zero. 


A number of methods for system modeling and parameter identification, in terms 
of effectiveness and complexity, are investigated in this thesis. To improve the system’s 
overall performance and stability under various operating conditions and faults, the 
results of the parameter identification are applied to the design of an adaptive controller 
for a power grid connected WECS. A basic view of the overall system is illustrated in the 


following figure. 









J 





Rotor Voltages 


General view of the proposed WECS. 


Xvi 


Since there are two independent control channels, one for the active and one for 
the reactive power, and they both have different requirements, the effectiveness of using 
either two independent adaptive controllers or a combination of one adaptive and one 
classical PI controller is considered. One can reason that the reactive power has simple 
control requirements since it has to be kept to a constant value. The active power, on the 
other hand, has to follow a profile for maximum torque, which changes with wind 


velocity and the aerodynamics of the turbine. 


To test and assess the effectiveness of the proposed techniques and their ability to 
adapt to changes in dynamics, a number of different cases have been simulated, such as 


changes in internal resistances and wind gusts. 


It was found that DMRAC can be effectively applied in a WECS, which can 
significantly improve the performance of the overall system. To that end, a physical 


implementation to validate the results of this thesis is recommended. 
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I. INTRODUCTION 


A. BACKGROUND 


The possibility of a shortage in fossil fuels in the near future in conjunction with 
the continuously increasing demand for energy [1], [2] and environmental concerns about 
high emissions has contributed to a high interest in renewable energy (RE) worldwide. 
Among the various sources of RE, one of the fastest growing and most promising is wind 
power. New wind power capacity added globally during 2010 reached 39 GW, which is 
more than any other renewable technology and over three times higher than the 11.5 GW 
of wind power added worldwide in 2005. As a result, existing capacity increased more 
than 24% relative to 2009, with total global capacity close to 198 GW by the end of 2010 
[3]. The growth of the global wind power installations over the past 15 years is 


demonstrated in Figure 1. 





1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 


Figure 1. | Global wind power installations, 1996 to 2010. From [1]. 


In the U.S., the wind industry started to grow in the 1980s; by 1986, the installed 
capacity in California was close to 90% of the global installed wind power. It later hit a 
plateau during an electricity restructuring period in the 1990s and regained momentum in 


the past decade [4]. 


The spread of wind energy (WE) installations indicates the need to develop wind 
turbines (WT) with higher efficiency. This involves technology advances that reduce the 
initial cost of the WT or allow the capture of more energy during operation, which results 
in a cost-of-energy (COE) decrease [5]. One approach that can lead to more effective WT 
is the design of more sophisticated controllers. The problem of designing an adaptive 
controller for a variable speed horizontal axis WT (HAWT) using a doubly-fed induction 
generator (DFIG) is addressed in this thesis. 


B. THESIS OBJECTIVES 


The objective of this thesis is to apply various methods of system identification to 
a DFIG and investigate the design of an adaptive controller for a power grid connected 
wind energy conversion system (WECS) to improve the system’s overall performance 
and stability under various operating conditions and faults. To achieve those goals, 
attempts are made to control both mechanical torque and reactive electric power 


independently using the rotor voltage components as a control input. 
C. RELATED WORK 


Various researchers over the past several years have addressed the problem of 
efficient control of WECSs. The complexity of the overall system and the random 
conditions under which WTs operate have led scientists to investigate different control 
approaches to improve performance and stability under various fault conditions. A novel 
control method for smoothing the stator active or reactive power ripple components under 
unbalanced grid voltage is proposed in [6]. In [7], the design of a multi-input multi-output 
(MIMO) sliding mode (SM) nonlinear control strategy for a WECS is presented. The 
resulting robust control law guarantees finite time convergence, whereas smoothing 


discontinuities at the basis of sliding mode control achieves chattering reduction. The 


zi! 


paper concludes that this type of control is a particularly suitable option to deal with 
electronically controlled variable speed operating WECS. Furthermore, SM control has 
proven to be robust with respect to system parameter variations and external disturbances 
as well as poorly known operating environments [8], [9]. In [10], direct model reference 
adaptive control (DMRAC) was used to design an adaptive collective pitch controller for 
a HAWT. The objective of the adaptive pitch controller was to regulate generator speed 
for wind speeds higher than the rated value for the generator and to reject step 
disturbances. This objective was accomplished by collectively pitching the turbine 
blades. The results were compared in simulations with a classical proportional integrator 
(PI) collective pitch controller. In the simulations, the adaptive pitch controller showed 


better speed regulation for high wind speeds. 


A different approach of optimal adaptive control is presented in [11] for a 
permanent magnet synchronous generator driven by a pulse width modulation (PWM) 
inverter. The controller focused on optimization of the system by minimizing Joule losses 
in the inverter and generator. To achieve this, the paper proposed the operation of the 


system at an operating point with minimum phase current, considering fixed torque. 


Other, more classical approaches include the independent control of active and 
reactive power using back-to-back PWM voltage source converters in the rotor circuit. A 
vector control scheme on the supply side PWM ensures independent control of active and 
reactive power drawn from the supply, while vector control of the rotor side PWM results 
in a wide speed - range operation. The control is achieved using PI controllers to enable 
maximum speed tracking for maximum power extraction from the wind. This approach, 
presented in [12], presents two different controls, a current-mode control and a speed- 
mode control. A similar approach was presented in [13], where a Stator-Flux Oriented 
vector control approach is deployed for both stator and rotor side converters to provide 
independent control of active and reactive power using a decoupled design based on 


internal model control. 


D. APPROACH 


In the first part of this research, the analytical design of an online identification 
algorithm is presented, with the intent of constructing a Simulink model of the physical 
system. In the second part, an adaptive control scheme to compensate for dynamic 
variations of the generator is proposed and its effectiveness is evaluated by computer 
simulations. In the simulations, various fault conditions are implemented, and the results 


are discussed subsequently. 
E. THESIS OUTLINE 


This thesis is organized as follows: the concept of the basic theory of WTs and 
DFIGs is covered in Chapter II, and the model created for the system in the Simulink 
environment is described. The theory of online system identification including modeling 
and simulation results for different estimation algorithms, is covered in Chapter II. A 
basic introduction to adaptive control theory and the design and simulation of two 
different adaptive controllers, along with simulation results under different operational 
conditions of the system, are presented in Chapter IV. A comparison with simulation 
results from other control approaches is also presented. Finally, the conclusions based on 
the results from the previous chapters, along with some recommendations for future 
research in the area of adaptive control for WTs, are presented in Chapter V. The 
Appendix includes Matlab code for initialization of the simulation of the controllers as 
well as the basic model blocks in the Simulink environment.II. WIND TURBINES — 
DFIG 


A. WIND TURBINE THEORY 


Wind turbines are designed to capture a portion of the wind's kinetic energy and 
convert it to usable energy [14] in the form of electricity. This section presents a brief 


explanation of the basic operating principles of WTs. 


1. Energy in the Wind 


Wind speed is a random process that depends on several factors. The main 
parameter that describes the ability to extract power is the mean wind speed for the 
specific area of interest. The probability distribution function (PDF) of the mean wind 
speed V,,, over an area is modeled as a Weibull distribution [14] of the form: 


P(V, y= H(t) oe (0.0) 


m 
CVC ’ 


where the parameters k and c are called the shape and scale factor, respectively, and are 


both positive real numbers. 


Figure 2. A volume of air flowing with speed V. 


If one takes a volume of air A4x flowing with speed V as shown in Figure 2, the 
kinetic energy stored in the wind is 


1 


E=5mv’, (0.0) 


k 


where m is the mass of the air. 


If we replace the mass in terms of the air density p and the volume as m= pAAx,, 


Equation (2.2) can be written as 


pve 5 piAAx)v?, (0.0) 
and from Equation (2.3) one can determine the wind energy per unit volume AAx =1 as 
f= ; pv. (0.0) 


Since the power is defined as the derivative of the energy with respect to time, 


one can derive an expression for the power of the wind: 


1 2 
2 AE 7 Oy sa cet 
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= 5 PAV". (0.0) 


2. Wind Turbine Aerodynamics — Disc Actuator Model 


Turbine aerodynamics describes the forces developed on a wind turbine by an 
airflow passing through it. The approach used to derive a model for the wind turbine is 
actuator disc theory [14]. In this theory, the WT is regarded as an actuator disc, a device 
that is designed to extract part of the kinetic energy of the wind. This theory also gives a 
value for the maximum efficiency limit. If one assumes the disc is immersed in an airflow 
as shown in Figure 3, then the upstream wind speed V has to be greater than the 


downstream wind speed V.. 





Figure 3. Actuator disc. After [14]. 


Furthermore, for the stream tube just enclosing the disc, the upstream cross 
sectional area A. is smaller than the disc area Ap, which in turn is smaller than the 
downstream cross-sectional area A... This happens because, by definition, due to the law 
of mass conservation, the mass flow must remain constant within the stream tube [14]. 


The above statement implies: 
AM 
——=constant , (0.0) 
At 


where AM is the air mass going through a surface during time At. Now one can derive the 
following expression: 


AM Ax 
—— = pA— = pAV =constant. 0.0 
At - At i“ ae 


Therefore, if one takes ideal points away from the turbine (— and + infinity) and 
on the turbine, one obtains: 


PA,.V = PA,V, = PA_..V.., (0.0) 


Since the air that passes through the disc undergoes a speed reduction equal to 
V—V.,,, this implies a reduction in its kinetic energy. This results in a force Fp developed 
by the actuator disc on the incident airflow, which is given by Newton's second law of 
motion as the product of mass M and acceleration a: 


VV) 
At 





F, = pA, Ax =(V-V_,) PAW, . (0.0) 


The force Fp originated by the pressure drop introduced by the actuator disc can 

also be expressed as: 
F, =(25 — Pp)Ap » (0.0) 
where, p; and p, are the values of the air pressure immediately before and after the disc, 


as shown in Figure 4. 


a 





Actuator disc 






Figure 4. _ Pressure around actuator disk. After [15]. 


Since Vp is smaller than V, one can write the speed at the disc as 
V,=(1-aV, (0.0) 
where one defines a, with 0<a@<1, which is referred to as the “axial flow interference 


factor.” Substituting Equation (2.11) in Equation (2.9), one gets the force Fp as: 


F, =(V—-V_,.)pA,V(1—a). (0.0) 


Equating Equations (2.10) and (2.12), one obtains: 
(V—V_,,)eV(1—@) =(p; —p). (0.0) 
For an airflow that is stationary, incompressible, and frictionless and with no 
external forces applied either upstream or downstream, the Bernoulli equation states that 
the total energy of the airflow remains constant. Given that the above conditions are 
satisfied, one can apply the Bernoulli equation from far upstream to just in front of the 


rotor and from just behind the rotor to far downstream [16]: 
io. “er <2 
(P, ~Po)= SPV" -V*.) . (0.0) 


Substituting Equation (2.14) in Equation (2.13), one finally gets 
V,=(1-2a)V, (0.0) 


and, since V_, has to be non-negative, Equation (2.15) gives us the upper bound of a, 


which is 0.5. 


Substituting Equation (2.15) in Equation (2.12), one gets an expression for the 


force that is applied from the airflow on the disc: 


F, =2pA,V’a(1-a). (0.0) 


The power that the actuator disc extracts from the wind is given by 
P, =F,V, =2pA,V°a(l—a)’. (0.0) 
This gives an idea of the maximum power available to the turbine from the flow 
of the wind. 


A conventional way of characterizing the ability of a wind turbine to capture wind 
energy is the power coefficient C,, [14], which is defined as the ratio of extracted power 


divided by the total wind power: 


P 
Ca, 0.0 
P (0.0) 


Substituting Equations (2.2) and (2.11) in (2.12), one gets the following 
expression for C,: 


_ 2pAV a(1-ay 


C 
0.5pA,V° 


=4a(1-a). (0.0) 





The power coefficient C, is a function of the axial flow interference factor a. The 
maximum theoretical value of the power coefficient can be found from Equation (2.19) 
by taking its first and second derivatives. The result, which is called the Betz limit, is 


C = 0.593 for a equal to 1/3 [17]. This limit is valid for all kinds of wind turbines, but 


p_max 
the actual efficiency is usually lower, and for most commercial wind turbines, a good 


value is between 0.40-0.50. 


In practice three effects that lead to a decrease in the maximum achievable 


power coefficient [18]: 
e Rotation of the wake behind the rotor. 
e Finite number of blades and associated tip losses. 


e Non zero aerodynamic drag. 


a Torque and Power on the Rotor 


The relation between the energy of the wind flowing through the windmill and the 
mechanical power available at the rotor is discussed in this section. The parameter B is 
the pitch angle of the turbine blades, and A the tip speed to wind speed ratio : 
28 

V 


A (0.0) 


A graphical representation of the parameters involved in Equation (2.20) is shown 


in Figure 5. 
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Figure 5. Turbine disc. 


Substituting the expression for the power coefficient Cp from Equation (2.19) into 


Equation (2.17), one gets the following expression for the aerodynamic wind power: 


C,(4,B) 1 
vP 


P, =2pA,V° 5 PAC, (A, BV? . (0.0) 





Substituting the expression for the disc area A, = R* into Equation (2.21), one 
gets: 
P, => parc, 2, BV" (0.0) 
From Equation (2.22), one can easily derive an expression for the aerodynamic 
torque produced by a turbine of radius R rotating at an angular rate Q, with wind speed V: 
P11 paR’C, (A, BV" 


© Oo 2 Q 


r r 


(0.0) 


In Equation (2.23), the substitution of Q, by its value extracted from Equation 


(2.20) provides the following expression: 
= 1 3 2 
= sree C,(2, BV (0.0) 


where the Cg is named the torque coefficient and is related to C, by 
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caf. (0.0) 


In Equations (2.21) and (2.22), it can be seen that Co and C, are written as 
functions of 2 and f. However, it ia assumed in this thesis that the pitch angle f=0, and 
the above coefficients are used as functions of the ratio A only. 

A parametric approximation of the torque coefficient Cg is given by the following 


expression in terms of 4 [7] 


Co(A) = o(F ie, (0.0) 


where a, b, c are dimensionless parameters that characterize the rotor of the wind turbine. 


If one substitutes Equation (2.26) in Equation (2.25), then 
Cp(A) =a(3-1Jee*. (0.0) 


The plot of a typical power coefficient curve versus / is shown in Figure 6. 
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Figure 6. Typical C, versus speed ratio J curve. 


The speed ratio A is a very important parameter in the dynamic model of the WT. 


Since one wants to operate as closely as possible to maximum efficiency to have 
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maximum power capture and given that the wind speed cannot be controlled and the 
radius of the blades is fixed, one needs to vary the rotational speed of the wind turbine to 
keep 4 close to its optimal value ,,. This is one of the advantages of variable speed WTs 
over fixed speed WTs, since they operate close to optimal efficiency only for a small 
range of wind speed values. 

For a typical wind speed profile with a Weibull distribution having a shape 
parameter 4 = 2 and scale parameter c = 8.5, the variable-speed turbine captures 2.3% 


more energy per year than the constant-speed turbine. This is considered a significant 


difference in the wind industry [18]. 


4. Regions of Control 


There are three different operational regions for the variable speed WTs, which 


are shown in Figure 7. 
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Figure 7. Example power curves for a typical commercial WT. From [5]. 


13 


In region 1, the section of the solid curve, usually located to the left of 4 to 5 m/s 
wind speed, the power available in the wind is low compared to the losses in the turbine 
system; therefore, the turbine is usually not operating in this region. The use of modern 


control strategies is not usually critical in region 1 [5]. 


Region 2, illustrated by the cubic section of the red curve, is generally for winds 
speeds between 4 to 5 m/s up to the WT's rated wind speed, usually between 12 to 15 
m/s. In this area, the WT is in operational mode, and it is desirable for it to capture as 
much power from the wind as possible. The goal in this region is to approach the Betz 


Limit curve as closely as possible. 


Region 3 is for operation above rated wind speed of the WT, i.e., the wind speed 
above which maximum peak power is produced. The goal of control strategies applied in 
this region is to limit the wind power captured so as not to exceed the designed electrical 


and mechanical load limits of the system. 


All research in this thesis centers on Region 2. 
B. DOUBLY-FED INDUCTION GENERATOR (DFIG) 


Induction machines are used in a wide variety of applications and are, without 


doubt, the workhorse of the electric power industry [19]. 


An induction motor is an electric motor that is driven by alternating current (AC). 
It consists of two basic parts: the stator, which is on the outside and has windings 
supplied with AC to produce a rotating magnetic field, and the rotor, which is attached to 
the output shaft on the inside. When balanced multi-phase currents flow through the 
stator windings, a rotating air-gap magnetomotive force (MMF) is produced, and, if the 
rotor rotates at different speed than the speed of the rotating MMF, balanced multi-phase 
currents are induced in the rotor windings. There are two different types of rotor 
windings: wound and squirrel cage. In Figure 8 the cross section of a typical three-phase 


induction machine with a squirrel cage rotor is illustrated. 
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Figure 8. | Induction machine cross section. From [20]. 


In single-fed induction generators, the rotor windings are short-circuited, while 
DFIGs take full advantage of both stator and rotor excitations. The rotor currents and 


voltages of the DFIG are at the basis of the control techniques developed in this thesis. 
1. Reference Frames — Dynamic Models 


Three phase signals (currents or voltages) can be represented either by the 
respective time domain components (called “abc”) or, more conveniently, by the 


expansion with respect to a rotating reference frame (called “qd0”). 


The equations converting time domain "abc" voltages into "qd0" components and 


vice versa are very standard and can be found in a number of sources [12], [19]. 


In the "qdO" reference frame for a symmetrical DFIG, we have the following 


equations for the flux linkages from [19]: 
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qd0Os >S lr, kaos 


PAs aor =~ =, odor + j(@- QO, )\A; wont 


—jaoda 


tigos © Vados 


(0.0) 
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where y, i, A refer to voltages, current, and flux linkages in the qd0 reference frame of the 
stator (subscript s) and the rotor (subscript 7). The parameter w is the electrical angular 
rate of the reference frame, determined by the line frequency (say 60 Hz in the USA); @, 
is the electrical frequency of the rotor circuit; 7; ,7, are the stator and rotor resistances; 


and p stands for the differential operator d/dt. 


Any qd0 variable, say f, can be expressed in a complex form as: 


ao =F, — if, +9. (0.0) 


Therefore, the flux linkages are defined as: 


Cine =[4,, —A 0] 


qs 


: ; : y (0.0) 
Agar) = [Au See 0] 
where the flux components are related to the currents by standard equations [19] 
Ag = laghgs + hay lige + hor) 
Aas =Lighgg + Li (i,, a ia) 
S a Li Ss 
ede Goas ak (0.0) 
Aer =i, wat ee at | 
her =Li hy a Li (ig, + lie ) 
Aor = Liles 


Expanding Equations (2.28) and combining them with Equations (2.31), one gets 
the following set of equations to characterize the DFIG [19]: 
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X), = Oki, (0.0) 
Xi = Ob, 
Xy =O Ly 


express all the impedances. In the above equations, wy, is the base or rated frequency (say 
the standard 60 Hz); @,; is the electrical rotational frequency of the rotor; and w- is the 
electrical rotational frequency of the stator. The variables Z, and L,, are the leakage and 
magnetizing inductances of the stator windings, and L’, is the leakage inductance for the 


rotor windings. Also, as mentioned, 7, and 7',.are the stator and the rotor resistances. 


The electromagnetic torque of the generator is given by the following expressions 


Es (2 )EE Jewtt isha)» (0.0) 


is i (3\(E (Aardar = Aariar) 
EAN (0.0) 


3\(P 74 
is (3) (Aasias ~ Aashis) 


In the above expressions, P is the number of magnetic poles of the generator, and 


[19], which are all equivalent: 


T. is positive for motor action and negative for generator action. 
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The total real and reactive power of a three-phase balanced DFIG in steady state 
conditions is given by [19]: 
3 
P, = 5 Mesto + Vilas) 


: . (0.0) 
Q, = 3 Mess ~ Vilis) 


2. DFIG Model 


Equations (2.28) clearly define a nonlinear dynamic system where the fluxes 4 are 
the states and the voltages are the control signals. Equivalently, from [19], one can 


redefine the states in terms of the gd0 current components of the stator and the rotor. 


The whole system is a nonlinear dynamic system and was implemented by 
Professor A. Julian in Simulink [21], as shown in Figure 9. The time domain voltage 
components vz, vs, Vc of the stator and the rotor are explicitly shown, together with the 
state variables ig, is, ic of stator and rotor, in the time domain. In addition, other physical 
quantities, like the torque on the load 7;, are shown, along with all relevant power 
components P, (active) and Q, (reactive). The model of the DFIG is described in more 


detail in [21]. 


Stator Voltages qdO 





Rotor Voltages abc 


DFIG_abc 


Figure 9. General view of the DFIG model. 
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This model is used extensively in the next chapters, both as a benchmark for 
simulating the dynamics of the DFIG and as the basis of the analytical design of the 


proposed algorithms. 
C. CHAPTER SUMMARY 


The basic theory behind WTs and DFIGs was presented in this chapter, and a 
description of the simulation model for a DFIG was provided. The quality of the model is 
of great importance, because the effectiveness of the controller depends on it. In the next 


chapter, three different algorithms of online parameter identification will be presented. 


19 


THIS PAGE INTENTIONALLY LEFT BLANK 


20 


Hil. ONLINE PARAMETER ESTIMATION 


A. INTRODUCTION 


The mathematical model described in the Chapter IJ, with its Simulink 
implementation, is representative of the DFIG physical system. In general, the various 
parameters defined (resistances and inductances) have to be estimated based on voltage 
and the current measurements of the physical system itself. The mathematical model can 


be used for control design and validation by computer simulation. 


The purpose of this chapter is to describe the methods used to achieve on-line 
estimation of the DFIG parameters. The essential idea behind on-line estimation is the 
comparison of the observed system response y(¢) with the output of a parameterized 
model whose structure is the same as that of the plant model [22]. In this case, the output 
vector y(t) consists of all the measured currents in the qd0 reference frame under different 


conditions of excitation through the input voltages. 


The parameter vector @(¢) contains estimates of all the parameters in the state 
space model in Chapter II, and it is adjusted continuously so __ that 


y(0,t) approaches y(t) as t increases. Under certain input conditions, y(0,t) being close 


to y(t) implies that Olt) is close to the unknown parameter vector 0(t) of the plant model. 


The online estimation procedure, therefore, involves three steps [22]: 
e Create an appropriate parameterization of the plant model; 


e Design an algorithm for generating and updating the parameter vector 


estimate; and 


e Properly excite the plant such that the parameter estimates approach the 


actual expected values as t>o. 


In general, the plant model of a system to be studied for identification should be 


expressed in the following form: 
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y(t) =" (t)0(t) (0.0) 


which can be explicitly written as: 
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where the nx1 vector y(t) consists of the observed variables; the mx1 vector @(¢) consists 
of the unknown parameters to be determined; and the nxm matrix ®'(t) consists of known 
functions that describe the plant model, which are called regression variables[23]. In this 


thesis, as one can see later on, n=6 and m=S. 


It is shown in [22], [23] that under conditions of persistency of excitation, the 
estimated parameters converge to the true parameters. Loosely speaking, this is satisfied 
when the signal is sufficiently excited so that the input-output behavior can be modeled 


by a unique parameter vector. 
B. SYSTEM EQUATIONS 


Based on what was discussed in the previous section, one needs to rewrite the 
dynamic equations in terms of parameters and regression vectors. In particular, the five 


quantities that form the parameter vector 0 to be estimated are the following: 


Bali WN, Ee deals (0.0) 


Ss is 


The vector y(t) of observations is defined by the stator and rotor voltages: 
T 
y(t) =v(t) = [Vac Wie. Vee Wipe Wi Vor | . (0.0) 


Substituting Equations (3.3) and (3.4) in Equation (2.32), one obtains: 
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where the matrix on the right of Equation (3.5) becomes ®’ (t) in Equation(3.2): 


v(t) =" (t)O(t), (0.0) 


and q@ is the angular velocity of the arbitrary reference frame that was chosen. 


C. ESTIMATORS 


In order to estimate the machine parameters vector O(t), three different 


identification approaches were used: 


e least mean squares (LMS), 
e recursive least Squares (RLS), and 
e recursive least squares with exponential forgetting (RLS-EF). 


The next section describes the three different estimators. 


1. Least Mean Squares (LMS) 


The first estimation algorithm is based on the Least Mean Squares (LMS) method. 
This was initially formulated at the end of the 18th century by the German mathematician 
and scientist K. F. Gauss to calculate the orbits of planets and asteroids. From [23], using 
the above method, the parameters of the model should be chosen in such a way that the 
sum of squares of the differences between the observed and computed values, multiplied 


by numbers that measure the degree of precision, is a minimum. 
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a. Equations Setup 


Based on the model described in Equation (3.1), one can define the 


prediction or estimation error e(t) [23], as the difference between y(t) and y(t): 
elt) =y(t)—yit) = y(t)- 0" (r)0(t). (0.0) 
In addition, one can define the parameter error A(t) as follows: 
Alt) = Alt) - O(t) (0.0) 


Thus, from Equations (3.1), (3.7) and (3.8), one can derive the following 


expression for e(t): 
eft)=0"(t)Az) (0.0) 


From [24], the general form of an identification algorithm is defined by a 


differential equation, called the update law, which is of the form: 


O(t) = F{y(t), e(t), A(t), (t)}. (0.0) 


By the standard steepest descent algorithm, one can obtain the following 


expression for the derivative of A(t) with respect to time: 
Olt) =-p@(t)e(t), (0.0) 


where 4 is an arbitrarily chosen strictly positive parameter called the adaptation gain and 
allows for a variable rate of adaptation of the parameters. The greater the value of w, the 
faster the convergence, but the response is more sensitive to noise and might induce 


oscillations in the system. 


To analyze the behavior of this system, one can define a Lyapunov 


function: 


v(Alt)) = = [Atel (0.0) 
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According to [25], for A(t) to converge, the Lyapunov function described 
in Equation (3.12) has to be positive definite, and its total derivative with respect to time 
must be negative semi-definite. By definition, it is positive definite, since in this case the 
parameters are positive real numbers and this case takes the norm of the parameter vector. 
The derivative with respect to time of V(0(t)) is given by: 


dV(O(t)) 
t 


_ 6" (t) alt). (0.0) 


V(A(t)) = 





Substituting the expression of Olt) from Equation (3.11) into Equation 
(3.13), one can obtain: 
V(O,t)= —10 (t)P(t)e(t) . (0.0) 
Since the prediction error e(t)is related to the parameter error by 
Equation (3.9), one can obtain: 
V(0,t)=—w (t)(t)®" (t)O(t) =—z| © (t)O(t) P<. (0.0) 
The outcome of Equation (3.15) is that V(0,t) is negative semi-definite, 


and as a result, 








A(t)| decreases monotonically with time. Under conditions of persistency 


of excitation, this implies that the parameter error actually converges to zero as: 
A(t)-—A(t) > 0. (0.0) 
Since the system parameter vector @(t)is assumed constant or slowly 


changing with time, one can assume that @(t)*0. Therefore, by differentiating Equation 


(3.8) and substituting Equation (3.11), one gets the following expression for the 


derivative of the parameter estimates: 
At) =—Olt) = HO (the(t). (0.0) 
Integrating Equation (3.17), one gets the expression needed for the 


estimates: 
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Alt) = (0) + | 4P(r)ly(z)— (2) O(z) Hz , (0.0) 


where the initial condition can be arbitrarily chosen based on any a priori knowledge of 
the plant. 


b. Implementation 
The simulation of this estimation method required the creation of a 


Simulink subsystem which was added to the model described in Chapter II for the DFIG. 


A general view of the estimator is shown in Figure 10. 







Q(t) Vector 


Estimator 


Parameters 


V(t) Vector 


Vit) 


Figure 10. General view of the LMS estimation model. 


The Simulink blocks that provide the vectors ®(t) and V(t) are common 
for all the different methods and were implemented using Equations (3.6). These are the 
equations for reference frame transformation described in [19]. The estimator block for 


the LMS is shown in more detail in Figure 11. 
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Figure 11. Detailed block diagram for the LMS estimator model. 
It can be seen in Figure 11 that instead of using the same adaptation gain yw 


for all the parameters, different values of « for the different channels were employed. The 


adaptation matrix in Figure 11 is a diagonal matrix of the form: 


HH, Oo O O O 
0 gw 0 OO 
u_Matrix=| 0 O 4, O O 
0 O0O O “w&, O 
o 0 OO O pel 








The reason behind this is to compensate for different scaling factors and 
bring all the signals to a similar order of magnitude. This seems to improve the 


convergence rate and prevent instability. 


The values of the non-zero terms in the “ matrix were chosen to create a 
relatively fast convergence without affecting the stability of the system. Large values of 
the parameters ju; lead to fast convergence, but the estimates are oscillating. Very small 
values lead to very slow but stable convergence. To illustrate, the values chosen for the 


above parameters that gave a stable output are shown in Table 1. 


va | 


Table 1. | Adaptation gains for the LMS simulation. 


























Hii Value 

Hi 130 

H22 820 

L133 0.007 

[44 0.38 

[s5 0.0135 
2. Recursive Least Squares (RLS) 


The second method was recursive least squares (RLS). 
a. Equations Setup 


This method is implemented in discrete time and recursively computes the 

parameter estimates A(n) by minimizing the quadratic error [23]: 
1d a 
V(0,n)= 5 WT.) —@'(iT,)Oy, (0.0) 
i=1 

with 7, as the sampling interval. 

From [26], the recursive estimate can be interpreted as a Kalman Filter for 
the process: 

O(n +1) = O(n) 


. ; (0.0) 
y(n) = ®"(n)6(n) + e(n) 


where y(n) @®(n) are given by Equation (3.4) and Equation (3.6), respectively, by 


substituting ¢ with n. 


The recursive equations for the RLS can then be derived from the Kalman 


filter equations used in [27]: 
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O(n +1) = O(n) +K(n)[y(n +1) —®" (n +1)A(n)] 
K(n) =P(n)®(n+1)[O' (n +1)P(n)\@(n +1) +! 
P(n +1) =P(n)—K(n)[®" (n +1)P(n)®(n +1) +11 *K’ (n) 


; (0.0) 


where the matrix @(n) has to be full rank, which means that it has to be nonsingular for 


all 7 >0. The initial condition for Equations (3.21) is: 


P(0)=0°L xm (0.0) 


where a is a positive number, much higher than the parameters, and m is the dimension of 


y(t), which in this case is six. 
b. Implementation 


The general outlook of the RLS estimator is shown in Figure 12. 






Parameters 


TIME DELAY 


Figure 12. General view of RLS estimator model. 


A more detailed view of the inside of the RLS block is shown in Figure 


13. The inputs to the block are @(n) and v(n), and the output is the parameter vector 


estimates O(n +1). 
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Figure 13. Detailed view of RLS estimator model. 


ae RLS With Exponential Forgetting (RLS-EF) 


In the above estimations, the parameters are assumed constant or slowly changing. 
Now consider the case of varying parameters from either environmental or operational 
conditions or a change in the parameters due to a fault in the system. The case that is of 
more interest is when the parameters are slowly changing in time. In this case, one needs 
to design an algorithm with a “forgetting factor” which discounts data exponentially 


further back in time. 
a. Equations Setup 


According to [23], the RLS-EF algorithm is obtained by changing the cost 
function of Equation (3.19) from the one that was used for the RLS algorithm to the 


following: 
ViOxt)=5 2) -O" (ONT (0.0) 
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where A is a parameter such that 0</ <1 and is called the forgetting factor. The meaning 
of Equation (3.23) is that one weights the last data by 1, with a factor of 2" for data that 
are n time units old. This means that one reduces exponentially the weight of each set of 
measurements. Thus, the equation for the RLS with exponential forgetting becomes: 

Olt +1) = Olt) + K(t)[y(t +1) — 0" (t + 1)(t)] 
= P(t) 


K()= Oe +O +a) oe+ ye (0.0) 
P(t+1)= ro. K(t)[®" (t +1) Mo ea KT (t) 
b. Implementation 


The general outlook of the RLS-EF estimator is similar to that shown in 


Figure 12. 


D. SIMULATION OF THE ESTIMATORS 


i; Equipment 


The three estimators described in the previous pages were simulated using 
Simulink software. A Matlab code was created to initialize the models. The code used for 
initialization of the simulation models in Simulink is presented in the Appendix. To 
determine how well the estimators perform, one has to choose a DFIG with known 
parameters for comparison. The DFIG that was used is the Lab-Volt model 8231, which 
is a 175 W, 120 V, 60 Hz, 4-pole machine. The system parameters were used to 
effectively model the system in Simulink. Both the stator and rotor of that machine 


consist of 3-phase wye-connected windings and a stator to rotor turns ratio N,/N,. = 10. 


To have an input to determine the values of the estimated parameters required use 
of the values that were measured in [28] after a DC test, a no load test, and a blocked 


rotor test on the generator mentioned above. The results are given in Table 2 
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Table 2. | Measured DFIG parameters. From [28]. 
































Parameters 0(t) Values 
Ee 12[Q] 
rp 15[Q] 
L, 0.0241[H] 
L, 0.0241[H] 
Lu 0.3342[H] 
2. Simulation Results 


After conducting a large number of simulations, varying the conditions and some 
parameters of the system, it turned out that only the first two methods gave satisfactory 
results. The simulation using RLS-EF method was leading the estimator to an unstable 


output after a period. The next section describes the results of the first two methods. 
a. Simulation Using LMS 


The results for the parameter convergence using LMS are shown in 
Figures 14, 15 and 16, where the initial condition for the estimator was set to zero. The 


sampling time for all the simulations was set equal to 0.1 ms. 
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Figure 14. 
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Figure 15. 
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Figure 16. Convergence of L;, and L;, using LMS from simulation. 


One can see that the parameter estimates converge to the expected values, 
but it takes a while for them to settle down. This is a result that depends on the non zero 
values of the adaptation matrix given in Table 1. If one sets the values of initial condition 


vector closer to the measured values, then the convergence is much faster. 
b. Simulation Using RLS 


The results of the simulation using RLS are shown in Figures 17, 18 and 


19; 
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Figure 19. Convergence of L;, and L;. using RLS from simulation. 


From Figures 17 to 19, one can conclude that the parameter convergence 
using RLS is almost instantaneous, but for the inductances, the estimates oscillate around 
the expected values. On the other hand, the estimates for the resistances are extremely 


close to the experimentally measured parameter values. 
E. CHAPTER SUMMARY 


In this chapter, different algorithms for system identification were presented 
analytically and simulated. Based on the results and given the fact that the nature of the 
system is relatively slow, one can say that the LMS is more suitable in this application 
since it performs well and provides relatively fast convergence. In the next chapter, a 


basic introduction to adaptive control theory and the design and simulation of two 
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different adaptive controller approaches, along with simulation results under different 
operational conditions of the system, will be presented. A comparison with simulation 


results obtained from classical control approaches will also be presented. 
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IV. ADAPTIVE CONTROL 


A. INTRODUCTION 


According to Webster's dictionary, the definition of adapt is "to change (oneself) 
so that one's behavior will conform to new or changed circumstances" [22]. According to 
another definition, an adaptive system is any physical system that has been designed with 
an adaptive viewpoint [23]. From the above expressions, one can say that an adaptive 
controller is a controller that has the ability to change its behavior to adapt to changes in 
the systems dynamics. That means that to be adaptive, the controller has to have 


adjustable parameters. 


Research on adaptive control systems was introduced during the 1950s, mostly for 
the design of high performance flight control systems. For many years, there has been 
skepticism among scientists about the use of adaptive control, since, in spite of the fact 
that the adaptation feature sounds very attractive, its performance and reliability may be 
questionable. In particular, when one applies it to a dynamic system, the designer should 
be aware that adaptation adds another loop, which might make the overall system 


unstable. Some important issues are the following [29]: 


e Because of internal modulation effects, fast adaptation produces a high 
frequency input. This input excites unmodeled dynamics easily and causes 
systematic instability. Typical examples from [30] have shown that under 
very mild conditions of unmodeled dynamics, an adaptive control system 


can easily become unstable. 


° There is no guarantee of the performance of the closed loop system during 


the transient response. 


In the last few years, new applications of adaptive control have started to appear, 


and some of these issues have been addressed with a degree of success [31]. 
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An adaptive control system can be thought as having two loops. One loop is a 
normal feedback loop, with the process and the controller. The second is the parameter 


adjustment loop. 


The system is modeled in terms of a vector of parameters @ that may be known or 


approximately known. The goal of the controller is to continuously produce a vector of 


estimates for the above parameters 6 and control the output vector. 
B. ADAPTIVE CONTROL THEORY 


There are many different approaches to adaptive control. The model reference 
adaptive control (MRAC) and, more specifically, the direct MRAC (DMRAC), is used in 
this thesis. In this control scheme, the system is modeled in terms of a vector of 
parameters 0. that may be known or approximately known. This parameter vector is 
updated by an adaptive law using input and output data. A general view of the direct 


MRAC is shown in Figure 20. 





Reference Model 


Wp) 






r(t) 


Controller 


C(@. ) 





z(t) 


Figure 20. Direct MRAC overview. After [31]. 


As shown in Figure 20, the output of the system z(t) is compared to the output of a 
reference model z,./(t). The reference model is chosen to generate the desired trajectory 
for the plant output to follow [22]. The controller parameters are updated continuously, 


depending on the output error: 


40 


z(t) a z(t) = _ ef (t) : (0.0) 
The goal of the controller is to continuously estimate the parameters giving a 
vector A(t) to drive the error of the output vector z(t)to zero. The rest of this chapter 


examines the application of this to the dynamic model of the DFIG. 
€. ADAPTIVE CONTROL OF DFIG 

Designing a control scheme required the evaluation of two different control 
approaches according to which part of the system should be adaptive: 


e First approach: adaptive control of torque and reactive power using as 


output the rotor currents i,, and ig,. 


e Second approach: adaptive control of torque and PI control of the reactive 


power using as output the rotor currents i,,and iz-. 
These approaches are described in the following sections and the results obtained 
from simulation are discussed. 


1. First Approach: Adaptive Control for Torque and Reactive Power 


In this approach, the controller was designed using adaptive control for both rotor 
currents components, i,, and ig, which form the z(t) vector. The overall design is 


described in the following section. 

a. Dynamic Model 

The voltage equations of a DFIG expressed in an arbitrary chosen qd0 
reference frame were described in Chapter II. 


For the purposes of this design, one should assume that the reference 
frame is determined by the line voltage of the stator. In this frame, the stator voltage is a 


constant with value: 


V,(C)=Vas(€) = Vat) = Jy (0.0) 


4] 


In addition, one can assume the stator resistance to be negligible; 1.e., 


r=0. (0.0) 


s 


Then, defining the rotational frequency of the frame @ to be the same as 


the electrical rotational frequency of the stator w,, one can write Equation (2.28) as: 


pa, =—jo, 


dqs 


IV, (0.0) 


The solution of Equation (4.4), if the initial condition is 1,(0)=“, isa 
a) 


constant given by 
A,(t) at. (0.0) 
a) 


Since the flux is complex /,(t) = 4,,(t)— j/,,(t), this yields 


Aglt)=—,constant (0.0) 


A,,(t)=0 


The fact that the flux is a constant in this reference frame leads to the 
definition of a constant magnetizing current i,,;. The goal of the rest of this section is to 
rewrite the dynamic model of the DFIG in terms of the rotor current components i,, and 


ig and the magnetizing current i,,; only. 


Therefore, using Equations (2.31) from Chapter II relating the fluxes with 


A.) [Flgtty ly di 
DleMhee or. ate Alls (0.0) 
H ly by thy [Li 


Substituting Equations (4.6) in (4.7), one gets: 


the currents, one can obtain: 
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V, 
(Gott etki, = 
(0.0) 


Sal 


(L,, +Ly digs + Li. 


Solving Equation (4.8b) for ig,, one gets: 





ae ee ee (0.0) 








: tes. tt VY, 
Liphms = Ags => bing = = bins =—* (0.0) 
- oly 


and substituting Equation (4.10) in Equation (4.8a) and solving for i,;, one obtains: 


V 
Liles — he hha bal 
- (0.0) 


M. 


i =< 
qs 
ee 


Bae 


is 





ms 


Now one can easily rewrite all the equations in terms of ins and i). 
Substituting Equations (4.9), (4.10) and (4.11) in Equations (4.7) and expanding the 


expressions for A,,andAg;, one obtains: 


Oi 6 ti (Ea | Mee 


qr 














_=L,(—*i,, Ll, +h). => 
q M L,, +L, be ‘i q I M"'q (0.0) 
hs Le 
Rot = = lig + (l= a. )(L ge Lal, qr 
ep (Ee ee | (ee 


and 


Near = Liphas + (L,, + Ly, Var > 





Ag = by | ia ee Ad |e (0.0) 
Is M 
L 2 
A, = (4 -—_—————__J(L, +L, )i 
dr ( (i, #4, it, 44,)" Ir ua Mar 
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The flux linkages equations can then be summarized from Equations (4.6), 
(4.12) and (4.13) as: 








V, 
Aas = a 
Ras =0 
ce ; (0.0) 
Qo = ing FO(L, + Ly Dit, 
qr a Fa ie ms ( Ir a) qr 
Rex = o(L,, + Lin dir 
where o is defined as a real parameter given by the following expression: 
ie 
7 <o<1. (0.0) 


,0 
(L,, o Ly \(L,, + Ly) 


In most cases, Ly is significantly larger than L;, and L;,, which means that 


o is very small. 


From Equations (2.28) in Chapter II, after ignoring the 0 components, 
which are negligible and defining 


O-0,=0, (0.0) 
one gets: 


pA, a lar 7 Oy, + Vor 


(0.0) 


Pag, = Wig, + Og, + Var 


r qr 


Substituting Equations (4.14) in Equations (4.17), one finally obtains the 


desired dynamic model of DFIG in terms of magnetizing and the rotor currents: 


: r, fae 
=H OE 
¥ ot, +t,)" ™ oft, +t,) 7 
; (0.0) 
; fF a ae Lv : 1 
Ply, = 


Se Op ee gee 
o(L, +L,,) o(L, +L,,)(L, +L,,) o(L, +L) 


Now, define a,a1,b the following parameters: 
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r 


aL Ly) 

= Li i 
o(L, +L,,)(L, +L)” 

<> a 

~ o(l, +L,,) 


(0.0) 





a, 


Substituting Equations (4.19) in Equations (4.18), one ends up with two 
simple first order models: 
pi,, =Agi_, — Diy, + bv, 


; na r. ‘ (0.0) 
Pig, =Agiy, + Oly, + OA, + bv, 


i V 
By defining the state vector Z Fr and the control vector v-| a and 
Vay Var 


choosing two arbitrary positive constants Gm, Bn, one can rewrite Equations (4.20) in 
terms of the state vector and a desired dynamic model: 


z(t) =p) at) "100 0] ! (0.0) 


with 
fos Qe 
ore-| : ‘| 
lar @ Vor 
(0.0) 
Ay + Qn 
0=| a, 
b- B,, 
The significance of Equation (4.22) is that by defining the control signal as 
follows 


G,,v{t)=—alt) lis | — D1 (t)O(t) +W, g(t) (0.0) 


qr 


45 


with A(t) an estimate of the system parameter vector O(t), the closed loop system 


becomes 
Z(t) =—a,, 2(t) + D" (t)O(t) + Wer (t) (0.0) 

where 
Alt)=0- Olt). (0.0) 


In the ideal case of perfect knowledge of the DFIG dynamics, the 
parameter error vector O(t) in Equation (4.25) is zero, and the state vector z(t) tracks the 


reference model: 


Zef (t) = —a,, ref (t) +E W oof (t) : (0.0) 
The tracking error 2(t) is defined as: 


z(t) =z(t)—Z,.(t). (0.0) 


Then, combining Equations (4.24) and (4.26) and substituting in Equation 


(4.27), one ends up with the following expression: 
Z(t) =—a,,z(t) +B" (t)O(t). (0.0) 


Now, estimating the controller parameters requires the development of a 


differential equation, called the update law. That differential equation is of the form: 


Olt) =Fty(t),elt), lt), O(t)}. (0.0) 


Using the steepest descent algorithm, we ensure that the parameter error 


decreases monotonically with time. In particular the definition of 6(t) is as follows: 


Alt) = pD(t)z(C) (0.0) 


where u>0O is the adaptation gain of the controller. The definition of the following 


Lyapunov function is: 
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V(6,t) = = [eel + 5 Pel. (0.0) 


As also stated in Chapter III, for the system to be stable, the Lyapunov 
function described in Equation (4.31) has to be positive definite, and its derivative with 
respect to time must be negative semi-definite. By definition, it is positive definite, since 
in this case the parameters are positive real numbers and the norm of the parameter 


vectors are used. The derivative with respect to time of V(O(t)) is 
V(A(t)) 220+ F000), (0.0) 
where, since O(t) 0, one can use the fact that: 
A(t) =—6(t). (0.0) 
In Equation (4.32), if one substitutes Equations (4.28) and (4.30), one gets: 


V(0,t)=2' (t)(-a,, 2(t) + ©" (t)0(t))—(W(¢)z(t) J"tt)=-a,, z(t). @.0) 


From Equation (4.34), V(0,t)is negative semi-definite for every value of 


time ¢. 


Equation (4.34) implies that: 




















V(t) <V(O)Vt 
Tie. Lie? ole, Ll iei2 > (0.0) 
=[2(t) +—|6(0)| <=||z(0) +—|6(0)| 
2 2u 2 2u 
which means, that if one chooses 
3(0) =z(0) => 2(0)=0, (0.0) 


And from Equation (4.35) 
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[ato 
a 


Equation (4.37) means that z(t) is bounded. 





ate)| Z (0.0) 


b. Simulation Model 


A Simulink model for the controller was created using the above equations 
for the system. The inputs to the controller are the following quantities, all expressed in 
the abc reference frame: stator voltages V,»,;, stator and rotor currents igocs and igney 
respectively, w, and the reference currents ij, ,.¢and ig ye The outputs of the controller 
are the rotor voltages that go into the DFIG model. A general block diagram of the 


controller is shown in Figure 21. 






PREPROCESS 


Parameter Estimation 





Reference Model 





Vaqdr_control 


CONTROL 


Figure 21. General view of the first controller in Simulink. 
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This model was created using Equations (4.18) to (4.37). A detailed view 
of the controller block, which gives the values of the control vector [Var Va | according 


to Equation (4.23), is shown in Figure 22. The output is then transformed back to the abc 


reference frame to supply the rotor with the proper control voltage. 





Vqdr_control 







omega_tild*[-idr iqr] 
Phi_trans(t) Matrix Vadr 


Multiply 
(Control Signal) 


Theta_hat(t) pi trans(t)*theta_hat(t) 


Figure 22. Control block for the first approach. 


The reference model shown in Figure 21 has a first order transfer function 


of the form 


Gn 





T(s)= (0.0) 


s+a,, 


The control goal is to independently control the rotor currents to minimize 
the reactive power of the DFIG and to produce the maximum real power for every 


different value of the wind speed respectively. 
From Equation (2.36), after substituting the values of V,; and Va; from 


Equation (4.2), one can obtain the following expression for the reactive power: 


(0.0) 


s qs? 


jay 
2 


which means that Q, is proportional to i,;, so to have Q, =0, ig; has to be equal to zero; 


L.€., 
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i 0. (0.0) 


qs_ref = 


From equation (4.11), substituting the above value for / from 


qs_ref 


Equation (4.40), one gets: 


tgp = Ves (0.0) 


which means that i,, should be equal to the magnetizing current, which is constant after 


the system reaches steady-state, and is given by Equation (4.10). 


On the other hand, keeping the power output maximum requires that the 
electric torque of the DFIG be equal to the aerodynamic torque that gives the maximum 


power output of the WT for each wind speed, which is given by [12] 


7. ref = Tie Satay = KG > (0.0) 
where 
1 slp (Aope) 
Kage = 5 PAR (0.0) 


and Q, is the rotational speed of the WT in revolutions per minute (RPM). This is given 
by 


1 2 


feet (0.0) 
Gear _Ratio P 


r re’ 


where Gear Ratio is the gear box ratio, P is the number of poles of the DFIG, and @, is 


the electrical rotational speed of the rotor in radians per second. 


The expression of the electric torque output of the DFIG is given by 
Equation (2.35). Substituting the values of the flux linkages from Equation (4.6), one gets 


3\(P Vie 3 
T, --(3)(F }ot% bls (0.0) 


and by substituting ig; from Equation (4.9) into Equation (4.45), one obtains 
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22 Oi» bok 


The goal is to track a reference torque, which requiresT, =, That 


e_ref* 
means that from Equations (4.44) and (4.46) one should set the reference current iy, that 


is driving the reference model to be 


ae (0.0) 


i = 
dr_ref (22) i” V, 
2/2) "L, +L, @, 


Running the model required the inclusion of a block that takes the current 








wind speed as input and gives the instant aerodynamic torque as output, which is the 


excitation of the DFIG. The WT turbine speed is then given by 

JO, = There —Tes (0.0) 
where J is the combined inertia of the WT and the DFIG rotor. In this thesis, the rotor 
inertia was neglected due to its very small value compared to the inertia of the WT. 


The comparison of the control approaches presented here with other 
controllers required the use of the parameters of the DFIG and the WT described in [7], 


which are shown in Table 3. 


Table 3. | WT and DFIG parameters. From [7]. 





























Wind Turbine DFIG 
R=3.8m Pnax=50 hp=37.285 kW 
J =3.362 kgm? P=4 
a=19.346 R=0.082 
b=9.4117 R=0.228 
c=20 L;=0.0355 Hy 
Gear_Ratio 16:1 L,=0.0355 Hy 
Aopt= 9.4 Ly=0.0347 Hy 
@p=_=2760 rad/s 
V; =380 V 
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Cc. Simulation Results 


After a large number of simulations, the system parameters that gave the 
best results are shown in Table 4. The sampling time for the simulation was chosen equal 
to 0.8 ms. For these values, the system was stable and tracked the reference model for 
wind speeds from 3.2 m/s to 10.8 m/s, which is a good range. This covered almost all the 
operational range of the WT. 














Table 4. Final simulation parameter values for first approach. 
Parameter Value 
din 100 
bn 250 
Mit 0.75 
b22 50 
[133 0.4 














In Table 4, it can be seen that the adaptation gain for the parameter 
estimation was chosen to be a diagonal matrix as seen in Chapter III to improve the 
convergence of the estimates. The initial condition for the estimates was chosen equal to 


70% of the expected values calculated using Equation (4.22). 


The wind speed profile used is shown in Figure 23 along with a plot of the 
absolute values of the actual output torque 7. and the reference torque T;. The wind 
speed profile was chosen in order to prove the performance of the controller even for 
large rapid changes in wind speed, greater than 5 m/s, which can be considered wind 


gusts according to [33]. 


The output torque is tracking the reference value for a range of wind 


speeds from 3.2 m/s to 10.7 m/s as shown in Figure 23. 
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Figure 23. Wind speed and torque during simulation of the first approach. 


Zooming in on the torque plot from time 90 s to 120 s in Figure 24, one 
can see that even for wind changes of 6.8 m/s, the difference between reference and 
actual output is less than 1.65 %. For lower wind speeds, the difference is always less 
than 0.7 %, which is a very satisfactory result. 
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Figure 24. Magnified torque plot during simulation of the first approach. 


The plots of the real power P. and reactive power Q, during the simulation are 


shown in Figures 25 and 26, respectively. 


a3 





Active Power(KW) 














10 20 30 40 50 60 70 80 90 100 110 120 
time (sec) 


Figure 25. Real power (P.) during simulation of the first approach. 
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Figure 26. Reactive power (Q-) during simulation of the first approach. 


In Figures 25 and 26, one can see that the ratio of QO. over Pe is very small 
and remains less than 0.2% for most the simulation time. This result indicates that the 
second goal of the design is also satisfied. The only time that Q, increases, but still with 
values less than 1% of P., is during extreme changes in wind speed, and that happens for 


a very short period. 


The rotor currents during the simulation are shown in Figure 27. 
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Figure 27. Comparison of the rotor currents i, and ig, with their reference values 
during simulation of the first approach. 


2. Second Approach: Adaptive Control of Torque, PI Control of the 
Reactive Power 


A second approach was attempted to simplify the controller with only one 
adaptive loop for torque control. The difference in this approach is that instead of 
controlling both i,; and ig, adaptively, a PI controller was used for i,; and the adaptive 


control for ig was kept. 
a. Dynamic Model of the Controller 


Equations (4.2) to (4.18) are valid for this approach. From Equation 


(4.11), solving for i,,, one gets 


Te ee ae (0.0) 


Substituting Equation (4.49) for the last one of Equations (4.18), 


one obtains 


a 


Pla, = 

















I wee by th x ae . 1 
=—-——_*+—_ jj, + Hi,,. — i,,)+@ Ae Vor 
o(L, +L) L o(L, +Ly,)(L,, +Ly,) o(L, +Ly) 
. (0.0) 
Oh! 3 ; 
fj getty 4 a a if Me Sag 
o(L, +L) c. Ch bbe) o(L,, + Ly) 
One defines the following parameters: 
r 
a, =-———— 
o(L, +L,,) 
a, abet 
Ly 
12 : (0.0) 
a, =(1 M Mis 
o(L, +L,,)(L,, + Ly) 
_ 1 
o(L, +L) 
Substituting Equations (4.51) in Equation (4.50), one gets 
Pig, =Agi4, —A,(@i,,) +a,0 + bv, . (0.0) 


Now Z,=/, and v=v, are called the state and control signal, 


respectively, which in this approach are scalars. Then, for any arbitrary constant an, Bm 


one can write Equation (4.52) in the form 


where 


and 


Z(t) =—a,, z(t) + B,v(t)+O"(t)0. (0.0) 





O'(t)=[iy —@ig O Vy, | (0.0) 
aA, +a, 
a, 
Bah (0.0) 
2 
b-£B,, 


In this case, the control signal is 
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Bnv(t)=—O"(t)O+W,(t). (0.0) 
Now substitute Equation (4.56) into Equation (4.53), and one gets: 
z(t) =—a,, z(t) +B" (t)0(t) +W,(t), (0.0) 
where A(t) is given by Equation (4.25). 
b. Simulation Model 


The general view of the controller in Simulink is shown in Figure 28. 
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PARAMETER ESTIMATION 
CONTROL 
th_tilde 


idr_ref 


inverse Ks transformation 
igs_ref 


REFERENCE MODEL 





[> 


Lz . 
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Figure 28. | General view of the second controller model. 


Equations (4.30) to (4.37), which were used in the first approach for the 


parameter estimates, were also used for this case. 
For igs, a PI controller was used with the transfer function 


Ty(s)=— =", (0.0) 


a7 


where Koy, K; are the gains of the PI controller. The input of this controller, the error 


between the reference and the actual measured value of i,; and the output is v¢-. 


Furthermore, Equations (4.42) to (4.48) that give the expression for ig, ref 


and Q, were again used here. 


In this approach, the output of the adaptive controller is vg, according to 
Equation (4.57); two low pass (LP) filters were added in the output of the adaptive part to 
filter higher frequency components that might affect the stability of the system. A 


detailed view of the control block diagram is shown in Figure 29. 








PID Controller 
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Transfer Fen 





igs_error 


Phi_trans Matrix 
Multiply 1/bm Transfer Fen 










Analog 
Filter Design 
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theta 
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Figure 29. Control block for the second approach. 


The cutoff frequency a, in the LP filter shown in Figure 29 was 100, and 
for the second filter, which is a 12" order Butterworth filter, the cutoff frequency was 10 


times higher. 
c. Simulation Results 


After a large number of simulations, the system parameters that gave the 


best results for this approach are those shown in Table 5. The sampling time for the 
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simulation was equal to 0.8 ms. For these values, the system was stable and tracked the 


reference model for wind speeds from 3.2 m/s to 10.8 m/s. 


Table 5. Final simulation parameter values for the second approach. 

















Parameter Value 
am 100 
bin 250 
Hi 100 











In Table 5, it can be seen that the adaptation gain for the parameter 
estimation in this case was a single value, which is about two orders of magnitude larger 
than the adaptation gains for the first approach given in Table 4. This value was chosen 
because comparable results could not be achieved between the two control approaches 
with smaller values of . The initial condition for the estimates was chosen equal to 70% 


of the expected values calculated using Equations (4.56). 


The wind speed profile is the same as the one used in the previous 
approach, which allows for a comparison between approaches in the next chapter. This 
plot, along with a plot of the absolute actual output torque 7, and the reference (optimal) 


torque T,.¢,is shown in Figure 30. 
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Figure 30. Wind speed and torque during simulation of the second approach. 
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From Figure 30, it can be seen that the output torque tracks the reference 
value pretty well for a range of wind speeds from 3.2 m/s to 10.7 m/s, but during the 
sudden wind changes, the actual output torque spikes for a few moments. A zoomed in 


view of the torque plot from time 97 s to 119 s is illustrated in Figure 31. 
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Figure 31. Magnified torque during simulation of the second approach. 


It can be seen in Figure 31 that for large and fast wind changes that 
simulate wind gusts, the torque tends to lose tracking for few moments. This is especially 
true when the wind decreases, e.g. at 100 s or at 117 s of simulation time, where the 


difference between the two values even reaches 100% but recovers quickly. 


The plots of the real power P, and reactive power Q, during the simulation 


of the second approach are shown in Figures 32 and 33, respectively. 
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Figure 32. Real power (P.) during simulation of the second approach. 
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Figure 33. Reactive power (Q.) during simulation of the second approach. 


It is can be seen in Figure 33 that Q, is quite high compared to P,. This is 


clearer if one zooms in on a part of the simulation, as shown in Figure 34. 
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Figure 34. Comparison between real and reactive power Pe and Qe, respectively, for 
part of the simulation during the second approach. 


Finally, the plots of the currents i,, and ig, are shown in Figure 35. 
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Figure 35. Comparison of the currents i,; and ig, with their reference values during 
simulation of the second approach. 


It can be seen in Figure 35 that iy, is tracking its reference value very well, 


but i,; seems to lose tracking during sudden wind speed changes. 
3s Performance Comparison 


To be able to reach a conclusion on the performance of the controllers designed in 
this study, a Simulink model of a classical approach for the control of a WECS was 
created. This is described in [12], again using the parameters of the system given in [7]. 


The model was created for [33]. 
Two different fault cases were simulated to compare the robustness of the control 
algorithms. These two are: 


e a wind gust, i.e., a sudden increase in the wind speed from 6 to 10 m/sec, 


° doubling of the value of the rotor resistance 7, to demonstrate what 
happens during a parameter variation. 
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The wind profile shown in Figure 36 was assumed. 
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Figure 36. Wind profile for comparison of the control approaches. 


The output torque of the three different control approaches during simulation is 


shown in Figure 37. 
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Figure 37. Comparison of output torque 7, and reference torque 7,.¢ during the 
simulation of a wind gust. 


Magnifying Figure 37 for about 10 s during the top of the wind gust provides 
Figure 38. 
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Figure 38. Magnified plot of output torque 7. and reference torque 7;.¢ during the 
simulation of a wind gust. 


From Figure 38, one can see that even though the performance of all the 
approaches considered is very good, the output torque of the first approach and the 
classical approach with the PI controllers are identical and almost identical, respectively, 
to the reference torque throughout the simulation, while the second approach shows some 
oscillation during sudden wind speed changes. Furthermore, the torque of the classical 
approach is slightly closer to the reference value, but has more overshoot compared to the 


first approach. 


The resistance of the rotor is a parameter that affects the shape of the torque 


output of the DFIG [19] and is illustrated in Figure 39. 
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igure 39. Typical DFIG torque versus time curve for different values of 7,. 


For the purposes of this thesis, in order to investigate the performance of the 


controllers that we proposed in previous sections, the rotor’s resistance 7, was increased 


to twice its initial value at t=15 s of the simulation, and all the control approaches where 


simulated. The output torque plots during the above change in r, from the simulations are 


shown in Figure 40. 
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Comparison of the variation of the output torque curves during simulation 
of a variation of r,. 
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As seen in Figure 40, 7. for the second approach has a high overshoot and even 
after three seconds does not completely track the reference torque 7, Between the other 
two cases, the first approach appears to have the best performance since it has a small 
spike in 7. but recovers faster than the PID approach. This can also be seen in Figure 41, 
where the speed of the DFIG for the three approaches is presented around the step change 
in r;; it is clear that the first adaptive approach gives the smoother speed curve among the 


three different cases. 
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Figure 41. Comparison of the variation of the DFIG speed during the simulation of a 
variation of 7. 


D. CHAPTER SUMMARY 


In this chapter, the analytical design, modeling and simulation of two MRAC 
approaches for a DFIG-based WECS were presented. The simulation results were 
compared with results from simulation of a more classical control approach for this kind 
of system, which was presented in [12]. The conclusions from these comparisons along 


with recommendations for future research are discussed in Chapter V. 
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Vv. CONCLUSIONS AND RECOMMENDATIONS 


A. SUMMARY - CONCLUSIONS 


The problem of the control of a doubly-fed induction generator based wind energy 
conversion system with fixed pitch horizontal axis wind turbine was addressed in this 
thesis. The adaptive control used was direct model reference adaptive control. This 
research effort addressed issues relative to the identification of the system’s parameters 
and the design of adaptive control techniques in order to compensate for uncertain or 


time-varying dynamics. 


In particular three system identification algorithms were presented based on 
standard techniques such as least mean squares (LMS) or recursive least squares (RLS) 
with and without exponential forgetting. The significance of these techniques is that a 
dynamic system of the DFIG can be properly parameterized on the basis of experimental 


results. 


Furthermore, two different adaptive control approaches were introduced and 
simulated for application in wind power generation. In this case it turns out that the 
parameterized dynamic model, based on the magnetic fluxes of stator and rotor, can be 
simplified considerably and only a few parameters have to be estimated. This is due to 
the stator being connected to the grid so that the stator voltage is constant in magnitude, 
frequency and phase. Given this characteristic, the magnetizing current becomes one of 
the parameters to be estimated, which results in considerable simplification of the 


dynamic model. 


These approaches have been tested in numerous simulations and seem to present 
an attractive solution to the problem of compensating for uncertain dynamics such as 


varying resistances and impedances in the circuitry. 
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B. RECOMMENDATIONS-FUTURE RESEARCH 


Extensive computer simulations demonstrate that the algorithms presented for 
adaptive control work as expected. Next step in this research will be the implementation 
of the proposed controllers in a laboratory environment using a direct current motor to 
simulate the wind turbine and building the controllers using Xilinx software in Simulink 
and field programmable gate arrays to establish bidirectional communication between the 


software and the equipment. 


This effort will involve not only actual real-time implementation but also the 
design of an appropriate model for the wind profile so that wind effects can be simulated 
by the generation of appropriate torque by the DC motor. The ultimate benefit of this 


research will be more efficient wind power generators. 
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APPENDIX: MATLAB CODES - SIMULINK MODEL 


In this Appendix, the code that was used for the initialization of the simulations, 


as well as the basic block of the Simulink models for the controllers are presented. 


A. MATLAB CODE FOR INITIALIZATION OF THE SIMULATIONS 


6% INITIAL CONDITIONS FILI 





PI 


Gaile; 
clear all; 





omega_b = 2*pi*60; 
tstep=.0005; 





tstop=120; 

twopiby3 = 2*pi/3; 
poles = 4; 

J=3.362; & kg m*2 
polesby2J3 = poles/2/J; 


GearRatio=16; 


SParameters from paper on nonlinear control 





rr0Q =0.228; % ohms 

rs=0.02; % ohms 

Xls=omega_b* (0.0355-0.0347); % ohms 

Xm =omega_b*0.0347; % ohms 

Xlr =omega_b* (0.0355-0.0347); % ohms 


Vs=380; %& Volts Grid phase voltage 
omega_e=omega_b; *%electrical frequency 











SParameters - Equivalent circuit 


is=Xls/omega_b; 
Lm=Xm/omega_b; 
ic=Xlr/omega_b; 
rsbyXls = rs/Xls; 























rrbyXlr = rr0/Xlr; 

Xaq = 1/(1/Xm+1/X1lst+1/Xlr) ; 

Xad = Xaq; 

XaqbyXls = Xaq/Xl1s; 

XaqbyXlr = Xaq/Xlr; 

XadbyXls = Xad/Xls; 

XadbyXlr = Xad/Xlr; 

V_phase = 380*sqrt(2)/sqrt(3); %Peak value used in Simulink for source 
model 
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SSteady state analysis 
omegar_ic = omega_b*.8 
0.001; 


j*Xm; 
23 = rr0/slip + j*Xlr; 
Ias = V_phase/sqrt (2) / 























3; 
Sstart at synchronous speed 


(Zl + Z2*2Z3/(22+2Z3)); 


Tar = -Ias* (22*23/ (Z2+2Z3))/23; 

Te = 3*poles/2*Xm/omega_b*real (j*conj (las) *lar) ; 
Iqse = sqrt (2)*real(lIas); 

Idse = -sqrt(2)*imag(lIas); 

Iqre = sqrt (2)*real(lar); 

Idre = -sqrt(2)*imag(lar); 

psi_mge = Xm* (Iqse+Iqre) ; 

psi_mde = Xm* (Idse+Idre) ; 

psi_qse = Iqse*Xls psi_mge; 

psi_dse = Idse*Xls psi_mde; 

psi_qre = Iqre*Xlr psi_mge; 

psi_dre = Idre*Xlr psi_mde; 

psi_gqsic=psi_qse; 

psi_dsic=psi_dse; 

psi_gric=psi_qre; 

psi_dric=psi_dre; 

SWIND 

rho=1.225; % air density 

lam=0.1:.1:10; 

C=(19.346./lam) .*(9.4117./lam -1).*exp(-20./lam) ; 


% plot(lam, C) 
lam_opt=6.4; 


C_opt=C (64); 





lamda=0.8; 


SMRAC parameters 
am=100.0; % 
bm=250; 

& Filter frequency 
ac=500; 


K_opt=rho*pi* (3.8%°5) *C_opt/ (2*lam_opt%2) ; 


this works 


sigma=1-(Lm*2) / ( (Lst+Lm) * (Lr+Lm) ) ; 


b=1/ (sigma* (Lr+l 
a0=-rr/ (sigma* ( 
c=50; 


° 


2 


% true values 
theta (1)=a0+am; 





Lm) ) ; 
Lxrt+Lm) ); 





% controller parameters 
im0=Vs/ (omega_e*Lm) ; % 


2 


theta (2)=Lm/ (Ls+Lm) ; 
theta (3)=-im0*a0*Lm/ (Ls+Lm) ; 
theta(4)=b*Lm/ (Ls+Lm) - bm; 


estimated magnetizing current 
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theta(5)=(Ls+Lm) /Lm; 
theta (6)=(1+(Lm%*2) / (sigma* (Lst+Lm) * (Lr+Lm) ) ) *im0; 
theta (7) =b-bm; 








° 


% Estimation bounds 


rho_theta=50; % parameter error in "percent" 
d_theta=abs (theta) *rho_theta/100; 
theta_min=theta-d_theta; 
theta_max=thetatd_theta; 


Q 


S$ initial estimates 
theta0=thetatd_theta.* (2*rand(1,7)-1); 





B. SIMULINK BLOCKS 


1. First Adaptive Control Approach 


a. Overall System View 











Manual Switch 





WIND SPEED PROFILE 


Aerodynamic lambda 
Torque (opt = 6.4) 










fal 


b. Aerodynamic Torque Input 








lambda 
1/2*p*pi*R5 


Constant 


U2 


Cc. General View of the Controller 
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d. Preprocess 


Computation of M(t) and 


@(t)T 





Fs and Fr Theta_c 





omega tild1 


omega_tild*[-idr iqr] 


omega_tild*[-idr iqr]2 


e. Parameter Estimation 


a3 


f Control 






omega_tild*[-idr iqr] Vaqdr_control 


Phi_trans(t) 


Theta hat(t) Phi_trans(t/"theta_hat(t) — 


(Control Signal) 





2. Second Approach 


The above Simulink blocks, in sections a, b, c, e, are common for both 


approaches. Those which are different for the second approach follow. 





a. Preprocess 
; Cc 
iqr 
Scope 
* * 
theta 
(ak) 
- = 
Phi_trans(t) 
( 3) 
idr 
ia) 4 
omega_tild Phi(t) 
omegs_tilde 
Terminator 
2 
C5) 
vqdr C4) 


Z(t) 





Computation of M(t) 
and M(t)" 
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b. Control 


iqs_error 


Multiply 








Phi_trans(t)*theta_c_hat(t) 


bm*Vqdor 


PID Controller 


10*[0.1, 0.1] 
s 
Transfer Fen 


1/bm 


ik 











Transfer Fen 
Analog 
Filter Design 


LP FILTERING 
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